The role of collective motion in examples of coarsening and self-assembly 
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The simplest prescription for building a patterned structure from its constituents is to add par- 
ticles, one at a time, to an appropriate template. However, self-organizing molecular and colloidal 
systems in nature can evolve in much more hierarchical ways. Specifically, constituents (or clusters 
of constituents) may aggregate to form clusters (or clusters of clusters) that serve as building blocks 
for later stages of assembly. Here we evaluate the character and consequences of such collective 
motion in a set of prototypical assembly processes. We do so using computer simulations in which 
a system's capacity for hierarchical dynamics can be controlled systematically. By explicitly allow- 
ing or suppressing collective motion, we quantify its effects. We find that coarsening within a two 
dimensional attractive lattice gas (and an analogous off-lattice model in three dimensions) is natu- 
rally dominated by collective motion over a broad range of temperatures and densities. Under such 
circumstances, cluster mobility inhibits the development of uniform coexisting phases, especially 
when macroscopic segregation is strongly favored by thermodynamics. By contrast, the assembly 
of model viral capsids is not frustrated but is instead facilitated by collective moves, which promote 
the orderly binding of intermediates consisting of several monomers. 



INTRODUCTION: SELF-ASSEMBLY AND 
COLLECTIVE MOTION 



Self-assembly refers to the generation of patterns or 
aggregates through the interaction of autonomous com- 
ponents [lj. The supramolecular self-assembly of thermo- 
dynamically stable structures plays a central role in biol- 
ogy, notably in the hybridization of nucleic acid strands, 
the organization of lipids to form cell membranes, and the 
assembly of proteins to form casings for viruses and other 
genetic material [5] . Biological components can also self- 
assemble when removed from their natural environment: 
protein complexes called chaperonins, for example, form 
large-scale sheet-like [3] and string-like [4] structures in 
vitro. Self-assembly is widespread in non-biological con- 
texts, from the myriad patterns formed by soap and oil 
in water to the ribbon-like structures that assemble from 
cobalt nanoparticles in solution [5]. The formation of 
such structures inspires the design of materials with novel 
properties, the chief goal of nanotechnology. 

The success of self-assembly is determined by an in- 
terplay of thermodynamics and dynamics. Components 
equipped with strong and highly directional interactions 
may stabilize the thermodynamically preferred structure, 
but are not guaranteed to assemble spontaneously into 
such a structure. Overly strong attractions may pre- 
vent structural relaxation (by impairing unbinding [51 [7] 
or 'reversibility' |8J), resulting in malformed aggregates. 
Components that must bind via precise alignment with 
a neighbor may take a prohibitively long time to do so, 
resulting in slow structural growth. This competition 
between the requirements for structural stability and ki- 
netic accessibility in general restricts viable assembly to 
small regions of parameter space, an idea illustrated in 
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FIG. 1: The competitive nature of self-assembly: an il- 
lustration of the antagonism between the requirements for 
structural stability and kinetic accessibility. We consider a 
schematic system of ovoid nanoparticles (mimicking, for ex- 
ample, inorganic nanorods), equipped with pairwise end-to- 
end interactions of strength J (darker shading of rod ends 
denotes greater interaction strength) and specificity, or in- 
verse angular tolerance, s (inversely proportional to shaded 
area). Stable structures of the required symmetry in general 
require strong, specific interactions (top left). However, such 
requirements tend to frustrate the kinetics of assembly (top 
right): too strong an interaction and structures fail to relax as 
they grow; too specific an interaction and productive binding 
events are rare. The competition between kinetics and ther- 
modynamics dictates the region of viable assembly (bottom) . 



Figure [T] 

Identifying and controlling regimes of viable assem- 
bly pushes the envelope of current experimental capabil- 
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ities. In support of such work, computer simulation pro- 
vides a powerful means of understanding self-assembly 
and its potential for creating new materials. Simulation 
permits exhaustive trials of model systems at little ex- 
pense [3 [Till EH [13 EH]- It can reveal both the nature 
of inter-component forces that lead to thermodynami- 
cally stable structures, and the dynamics through which 
such components associate jSJ |5J [H] ■ When assessing the 
assembly properties of a given model it is desirable to 
evolve that model in order to approximate the dynamics 
that the corresponding physical system would execute. 
Molecular dynamics [TH [TS] algorithms evolve compo- 
nents according to Newton's laws of motion, and so are a 
natural choice for simulating particle systems. However, 
self-assembling components generally possess anisotropic 
interactions of maximal strength much greater than k^T 
and range much less than a particle diameter. Strong, 
short-ranged interactions place stringent limits on the 
maximum integration time step able to preserve numeri- 
cal stability. Under such conditions, simulations of large- 
scale assemblies are very time-consuming, forcing the 
simulator to choose between focusing on dynamics on 
smaller scales or starting from forcibly equilibrated sam- 
ples. 

One way to circumvent this problem is to use a coarse- 
grained dynamical procedure to move particles accord- 
ing to potential energy gradients without explicitly in- 
tegrating equations of motion. The Monte Carlo tech- 
nique provides a flexible framework in which to do 
so [HI HQ [T71 HH1 HH HO]. However, conventional Monte 
Carlo techniques involve sequential moves of individual 
particles and so neglect the correlated motion of parti- 
cles on timescales less than the fundamental discrete time 
step At (the time corresponding to a typical discrete par- 
ticle displacement) . For some systems this neglect of col- 
lective motion appears to be unimportant: in Ref. |20j . 
for example, molecular dynamics results were reproduced 
using a single-particle Monte Carlo protocol. For many 
systems, however, chiefly those whose constituents pos- 
sess interactions whose strengths vary strongly with angle 
or distance, neglecting motion correlated on timescales 
less than At leads to unphysical relaxation, particularly 
at long times. 

To address such problems, 'cluster' algorithms have 
been used extensively to effect correlated or collective 
motion [2I1[11[131[M1[2S1[2S1I1Z1I1H1I1S]-In general, such 
algorithms identify collections of particles to be moved 
in concert by recursively 'linking' particles according to 
a set of criteria, such as the pairwise energy or degree of 
proximity of neighboring particles. One such algorithm, 
the 'virtual-move' Monte Carlo (VMMC) procedure of 
Ref. [30] . is designed to effect correlated displacements 
and rotations according to potential energy gradients or 
forces experienced under 'virtual' moves of neighboring 
particles. We describe this idea in the following section 
and in the Appendix. This procedure reduces to single- 
particle motion when particles experience small energy 
changes on timescale At, but effects collective motion on 



arbitrarily large lengthscales when particles experience 
large changes in energy on this timescale. 

In this paper we use VMMC to study the qualita- 
tive effect of collective motion on self-assembly by explic- 
itly allowing or suppressing collective moves of particles. 
In Section |III| we consider the 2d attractive lattice gas, 
which coarsens upon a temperature quench via the self- 
assembly of the homogeneous phase. We find a range 
of temperatures at which inter-particle forces are large 
enough to encourage assembly but not so large that mo- 
tion is strongly correlated on timescales less than At. 
In this case assembly is driven by single-particle binding 
and unbinding events, and little qualitative effect is ob- 
served upon accounting explicitly for correlated motion. 
However, at low temperatures we find that motion on a 
timescale At is strongly correlated according to poten- 
tial energy gradients. Neglecting or allowing explicit col- 
lective motion under these conditions selects drastically 
different fates for the system: strongly collective motion 
impairs assembly via the formation of kinetic traps asso- 
ciated with the binding of large clusters. In the language 
of coarsening, these kinetic traps represent the arrest of 
phase separation by gelation. Suppressing collective mo- 
tion suppresses gelation. We present a simple argument 
designed to estimate the importance of collective motion 
in the space of temperature and particle concentration. 

In Section II VI we examine a three-dimensional off- 
lattice system of hard spheres with isotropic pairwise 
square- well interactions, a model of strongly-associating 
colloids. Here, as for the lattice gas, collective motion 
tends only to impair assembly (or promote gelation) by 
inducing awkward binding events between large clusters. 
However, in Section [VJ we study a model of viral cap- 
sid assembly that displays qualitatively different behav- 
ior. Inter-particle forces are sufficiently strong that sub- 
stantial correlated motion emerges on the fundamental 
timescale At, but the geometry of inter-particle associa- 
tion is such that these collective motions improve assem- 
bly by inducing productive collisions between small in- 
termediates larger than monomers. Suppressing this cor- 
related motion slows assembly, but does not strongly im- 
pact the final capsid yield at most thermodynamic states. 
Collective motion therefore appears to play a qualita- 
tively different role within different model systems: it 
drives the formation of kinetic traps and facilitates or- 
derly growth. 



II. A 'VIRTUAL-MOVE' MONTE CARLO 
ALGORITHM 

Here we summarize the virtual-move Monte Carlo al- 
gorithm. We consider a <i-dimensional collection of N 
particles equipped with pairwise interactions. This al- 
gorithm is a dynamic procedure designed to identify, on 
the basis of potential energy gradients explored on a fun- 
damental timescale At, the extent to which the motion 
of one particle is correlated with that of its neighbors. 
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FIG. 2: Illustration of the VMMC procedure applied to a col- 
lection of pairwise-interacting oval nanoparticles. (a) Starting 
in microstate fi we pick a seed particle (shaded) from which 
to 'grow' a pseudocluster, together with a virtual move map 
(denoted by the arrow). The pseudocluster is a set of linked 
particles that will experience a trial move. The virtual move 
map is both a device for computing neighboring potential 
energy gradients, and defines the move that the pseudoclus- 
ter will execute. The pseudocluster growth procedure is an 
iterative linking scheme akin to the Swendsen-Wang (SW) 
algorithm |21J : we propose links between the seed and all 
particles with which it interacts, and continue iteratively un- 
til links have been proposed between the pseudocluster and 
all particles with which its constituents interact. Our pro- 
cedure differs from the SW algorithm in that here links are 
formed not on the basis of pairwise energies, but instead on 
the basis of pairwise energy gradients. These gradients are 
computed by executing virtual moves of shaded particles (ex- 
ample shown in lower panel), and recording neighboring pair- 
wise energies before and after those moves. Conditioning the 
acceptance criterion upon reverse virtual moves (not shown) 
ensures that detailed balance is preserved. The linking crite- 
rion is similar to that of the geometric cluster algorithm of 
Liu and Luijten [26] , although here we propose links only be- 
tween particles that interact in the initial configuration. We 
seek to effect dynamically realistic local and collective motion, 
rather than the collective and nonlocal equilibration-speeding 
motion of Ref. [26] . With the pseudocluster so defined (b) we 
displace it according to the virtual move map, resulting in mi- 
crostate v (c) . We then evaluate the Monte Carlo acceptance 
criterion (see text) and accept or reject the move. 



If so correlated, this motion is effected with a frequency 
designed to approximate a physical dynamics. An illus- 
tration of the procedure is given in Figure [2] We start 
in microstate \i. We define a 'pseudocluster' C (a set of 
particles to be moved in concert) by choosing as its first 
member a 'seed' particle i. We link the seed to a neigh- 
bor j with probability Pij(u — > u), which in general de- 
pends on a 'virtual' move of i that defines a notional new 
microstate, v. Particles linked to members of the pseu- 
docluster join the pseudocluster. We proceed iteratively, 
until no more members are added to the pseudocluster. 



We accept the move fi — > v with probability 
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Here T>(C) < 1 is a factor we impose in order to modulate 
the diffusivity of pseudoclusters according to size; E a is 
the energy of the system in microstate a; qij = 1 — pij 
is the probability of not linking particles i and j; and 
1Z denotes a particular realization of formed and failed 
links. The link-forming procedure is aborted in situ if 
the pseudocluster size nc exceeds a specified cutoff n c , 
the smallest integer larger than Here £ is a random 
variable drawn uniformly from the interval [0,1]. The 
subsequent rejection of the move is enforced by the fac- 
tor 9(ri c — nc). This rejection procedure ensures that all 
particles experience proposed moves with approximately 
equal frequency. The products over q variables quantify 
the probabilities of not forming proposed links, internal 
to and external to C (for forward and reverse moves), 
while the products over p variables quantify the proba- 
bilities of linking together members of the pseudocluster. 
We choose to link particles i and j with a probability 
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that depends on a virtual move (translation or rotation) 
of i relative to j. Here e{i,j) is the pairwise energy of the 
bond ij in microstate /i, and e(i',j) is the bond energy 
following the virtual move of i. The factor is unity if 
i and j interact in microstate /i, and zero otherwise; the 
factor O (n c — nc) terminates link formation if n c > nc- 
Linking particles in this fashion ensures that neighbors 
exert mutual forces proportional to the gradient of their 
pairwise energies, with motion 'linked' or correlated if 
(3e(i',j) — (3e(i,j) is large. Equation ^ implies that the 
acceptance rate ([!]) reduces to 
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The label {ij) n — identifies particle pairs that start (//) 
in a nonintcracting configuration and end (v) with pos- 
itive energy of interaction (overlapping), or which start 
(//) in an overlapping configuration and end (V) in a non- 
interacting one. The final line of Equation d3| accounts 
explicitly for links (ij)e and failed links {ij)t internal to 
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the pseudocluster (note that the q factors internal to the 
pseudocluster were omitted in error in Equations (11) 
and (13) of Ref. [3D]). In Appendix A we outline a pro- 
cedure in which this acceptance rate is simplified by 'sym- 
metrizing' link formation at the level of link generation. 



III. A PROTOTYPE OF SELF-ASSEMBLY: THE 
ATTRACTIVE LATTICE GAS 

In this section we consider the quench-driven coarsen- 
ing of the 2d attractive lattice gas with conserved particle 
number. Lattice gas models, whose thermodynamics can 
be related to the Ising model [3TJ [32], are used to cari- 
cature a diverse range of physical systems, from binary 
metallic alloys to solvent-mediated nanoparticle aggre- 
gation 33J . We regard coarsening within the attractive 
lattice gas as a prototype of self-assembly in which attrac- 
tive interactions drive the organization of a homogeneous 
phase. Our aim is to identify the range of temperatures 
and particle concentrations where collective motion (mo- 
tion correlated on a timescale less than the discrete time 
step At) strongly influences assembly. 

Many authors have studied the dynamics of the lat- 
tice gas using single-particle Monte Carlo algorithms, in- 
ducing transitions between microstates by moving a sin- 
gle particle to an unoccupied nearest-neighbor site. The 
assumption underlying these studies is that sequential 
moves of single particles represent a good approximation 
of the dynamics that the corresponding physical system 
might execute. This assumption is likely accurate for 
systems that relax via transport of mass from interfaces 
of high to low curvature by diffusion through the inter- 
vening medium [34]. Scaling arguments and simulations 
based on this physical picture predict the typical domain 
size L to grow as L(t) ~ t 1/3 [251134], 

However, in many settings, the transport of mass be- 
tween domains by the evaporation and diffusion of con- 
stituent monomers is not the only possible mode of re- 
laxation. Indeed, we frequently encounter the concerted 
motion of domains of one phase in another: witness the 
relaxation of polymers in solution or the the diffusion of 
nanocrystal aggregates on graphite [35]. To model such 
behavior within the lattice gas, we must explicitly ac- 
count for motion correlated on timescales of order At. 

We consider a collection of N particles on a simple 
square lattice of V — L 2 sites. Two particles may not 
occupy the same site, and interact with binding energy 
— eb when occupying nearest-neighbor sites. We express 
temperature T in units of et,/&B ■ Particles are dispersed 
randomly on the lattice with concentration O = N/V 
and are evolved using virtual-move Monte Carlo transla- 
tions. We enforce a scaling of the diffusion constant of 
■nT a for clusters of size n > 1. We focus on the difference 
between the case a = oo, corresponding to single-particle 
moves, and a = 1, corresponding to diffusion akin to 
Brownian dynamics. 

The simplicity of the model allows us to estimate the 



values of temperature T and particle concentration 0o 
for which we expect motion correlated on a timescale 
At to be important. We present this argument for the 
c?-dimcnsional hypercubic lattice gas. We assume for 
simplicity that growing clusters are compact (as a two- 
dimensional illustration, consider the left half of Figure [6] 
top row). We wish to estimate the timescale on which 
two clusters of size n encounter each other as a result 
of their collective diffusion, T eric (n, n), and compare this 
estimate with the timescale on which such clusters ex- 
change mass via evaporation (Ostwald 'ripening' [36 ), 
Tevap(^)- If T evap (n) exceeds T enc (n), we expect collective 
modes of motion to significantly influence assembly. 

Within the virtual-move algorithm the probability of 
whole-cluster motion is approximately 

P^oUn)=P^ Kn - 1) ^, (4) 

where pii n k = 1 — e~' 3Ae is the probability of linking two 
particles following a virtual move, and Ae = > is the 
change in energy resulting from separating those parti- 
cles. We impose a factor of n -1 , which ensures that par- 
ticles suffer attempted moves with approximately equal 
frequencies, and a factor of n~ a , to account for our cho- 
sen diffusion constant. The number £ reflects the ef- 
ficiency with which the recursive algorithm forms links 
within a cluster. This number is approximately l/(2d) 
for the version of the algorithm described in the main 
text, and is zero for the version of the algorithm described 
in the appendix. This distinction is unimportant at low 
temperature, and we shall for simplicity take £ = 0. 

The probability that a single particle breaks away from 
its host cluster is 

Punbind(n) « e~^« Ae g n , (5) 

where 

n (d-\)/d 

Qn = 2d (6) 

n 

is a geometric factor quantifying the likelihood that a 
chosen particle lies on the surface of a cluster (we expect 
this approximation to be reasonable for n > 10), and z n 
is the typical coordination number of a particle on the 
surface of a cluster of size n. 



In d= 2 we have that p who ic(n) cx e-C"- 1 )! ln »^| n -i-« ; 
while p U nbind( n ) ^ n" 1 ^ 2 ■ For sufficiently large n we 
therefore expect single-particle binding events to domi- 
nate. To determine the cluster size n* at which unbinding 
and collective motion are equally likely, we equate Equa- 
tions (g) and 

(n*-l)lnpi illk = -f3z n * Ae+ln(2d)+ (^-7^+0) Inn*. 
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FIG. 3: Estimate of the size n* of the largest cluster for which 
collective motion is at least as probable as the unbinding of its 
constituent monomers, as a function of particle binding en- 
ergy Ae. Estimate is derived from a simple argument (Equa- 
tion j7|) pertaining to the 2d lattice gas, assuming constant 
coordination number z n — 1,2 and 3. For the two larger val- 
ues of z„, n* is appreciable even for relatively modest values 
of Ae. 



Clusters smaller than n* will likely move as a whole. To 
estimate n* at low temperatures, we ignore the term log- 
arithmic in n* to obtain 



n*(T) 



Ae 



Ae/(fc B T) 



(8) 



In Figure [3] we plot n* obtained from numerical solution 
of Equation ^ for d — 2, a = 1, and three values of z n . 
We observe that n* is large even for relatively modest 
values of Ae/(/cBT) when z n > 2. 

We next estimate the timescale upon which two clus- 
ters of size n <n* collide. The cluster diffusion constant 
is 



D(n) = D Q np who i e (n), 



(9) 



where Dq is the diffusion constant of a monomer. The 
timescale upon which two clusters of size n encounter 
each other through diffusion is approximately 



: (n,n) = -D{n)- l i en {nf, 



where 



eff(n) = (<f> 1 



-l/rf 



2a rP 



(10) 



(11) 



is a measure of the distance separating clusters of size n. 
Here a = 1/2 is the monomer radius, and 7 is an exponent 
measuring the increase in distance between structures 
due to clustering (our simulations indicate that 7 f=a 0.5 
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FIG. 4: Estimate of the regime of collective motion for the 2d 
attractive lattice gas for the given values of cluster size n and 
diffusion exponent a. We show the intersection of the 'clus- 
todal' (locus of Vcol\{n) = 1) with surfaces of constant cluster 
size n as a function of T and <^o- For given n and a, collec- 
tive motion dominates the behavior of clusters of size n (or 
smaller) in the region below the corresponding line. Hence, 
collective motion dominates this kinetic phase diagram for a 
broad range of temperatures and particle concentrations, with 
quantitative differences evident upon changing the exponent 
a (governing cluster diffusivity) from 1/2 to 1. We assume 
constant coordination number z n = 2.5. 



in d = 2). For brevity we write Iq 



1 n a ~ 



(1 



2a. Then 
(12) 



We compare the cluster-cluster encounter timescale with 
the timescale T eV ap(^) required for monomers to unbind 
from clusters of size n and encounter other structures of 



size n: 



1 cvap 
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Tunbind(n-) 
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(13) 



At low temperature the unbinding timescale is much 
larger than the timescale for diffusion of a single particle 
between clusters. 



For a cluster of size n, Equations ( 12 1 and ([13| quantify 
the respective timescales for mass transport by collective 
motion, and for the unbinding and diffusion of monomers. 
We view the ratio V co ii ( n ) = 7"evap(^)/T e nc("-) as a mea- 
sure of the propensity for collective motion at a given 
scale n. The self-assembly of the homogeneous phase 
proceeds in stages via the appearance of structures of 
size n. At each stage, we expect collective motion to be 
important if P co ii(^) is greater than unity. At high tem- 
perature, single-particle motion dominates: V C oii(p) ^ 
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FIG. 5: Top row: Kinetic behavior of individual clusters 
within the 2d attractive lattice gas at different T. Starting 
from disordered configurations at concentration 0o = 0.1 we 
evolve each system according to virtual-move Monte Carlo 
dynamics with D(n) oc n . We 'capture' clusters of size 
Umax = 10 or 30 as they develop, and subject these clusters 
repeatedly to the same collective-move Monte Carlo proce- 
dure, recording (but not making) accepted moves. We plot 
the probability of the concerted motion of a sub-cluster of 
size 1 < n < n max , P m0 vc(w), normalizing data by setting 
fmovc(l) = 1. At the highest temperature the probability 
of monomer unbinding, P m0 ve(l), is rapid relative to that of 
whole-cluster motion, P m0 vc(w max ), allowing assembled struc- 
tures to relax as they grow. Note however that correlated 
motions contribute measurably to relaxation dynamics even 
at high temperature. At lower temperatures, collective mo- 
tion of the whole cluster predominates. Bottom row: typical 
30-member clusters obtained at the three temperatures con- 
sidered. 



« = 1 



a = oo 



equil. 



a - 1 



a = <x> 



equil. 



ft -o.i 




<fc=0-4 



FIG. 6: Configurations of the lattice gas at fixed time follow- 
ing a quench from a disordered state, evolved using Monte 
Carlo dynamics with cluster diffusivity D(n) oc n~ a and 
particle concentrations 4>o = 0.1 (top panel) and 0.4 (bot- 
tom panel). The equilibrium states (obtained using nonlo- 
cal moves, panel 'equil.') are in all cases phase-separated 
configurations. Top panel: At high temperature, collective 
(a — 1)- and single-particle (a = oo) motion give rise to vi- 
sually similar assembly behavior, while at low temperature 
collective motion drives the formation of gel-like kinetically 
trapped structures. Gelation is avoided by suppressing col- 
lective motion. Similar behavior is seen at higher concentra- 
tions (bottom panel), although collective motion at the lower 
temperatures is hindered sterically. Coarsening images cap- 
tured after 2 million MC steps per particle (1.5 million steps 
for T = 1/3 and <j> = 0.4). 



(/3Ae)" 1 — > when /3Ae — > 0. At low temperature, in- 
stead, collective motion dominates assembly up to large 
values of n: V co n(n) ~ n -w/d- a -2 le f3z n Ae £ - 2_ For the 

choices 7=1/2 and d — 2, we expect collective mo- 
tion to be a more effective means of mass transport than 
single-particle unbinding for clusters of size less than 

nmM ~ [4 2 exp (f3(z n ) Ae)] 2/(2 " +3) , (14) 

which can be very large at low temperature and low 
to moderate densities. Under these conditions, whole- 
cluster motion will dominate the system's assembly dy- 
namics on large length and timescales. 

The regimes of single-particle-dominated motion and 
cluster-dominated motion are separated by the locus of 
points (</>q,T*), defined for a given n by the equation 
T'coiiM = 1- We call this locus the 'clustodal', by anal- 
ogy with 'binodal' (locus of phase equilibrium), and 'spin- 
odal' (locus of the onset of spontaneous decomposition). 



In Figure [3] we show a plot of the propensity for collec- 
tive motion, V co ii(n), in the {4>q,T) plane, based on our 
simple argument. We consider two values of a (quan- 
tifying the cluster diffusion rate). The intersections of 
the clustodal with surfaces of constant n indicate where 
in state space collective motion dominates assembly dy- 
namics on a scale n. Thus for a thermodynamic state 
(4>o = 0.1, T = 0.2), say, our estimate suggests that as- 
sembly dynamics (for a = 1) for clusters of size 30 (but 
not size 300) is dominated by motion correlated on times 
less than the discrete time step At. 

In simulations of the lattice gas we indeed find that col- 
lective motion is important for a broad range of temper- 
atures and particle concentrations. In Figure[5]we offer a 
microscopic perspective on the meaning of the clustodal 
by analyzing the (averaged) motion of individual lattice 
gas clusters of sizes n max = 10 and 30 at three temper- 
atures. Starting from randomly dispersed monomers at 
concentration 4>q — 0.1, we capture clusters of size n max 
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FIG. 7: Growth of domain size t with Monte Carlo time t 
within the attractive lattice gas, for cluster diffusivity a = 1 
or a = oo at particle concentration ifo = 0.1 on a lattice of 
size V = 256 2 . We perform quenches to T = 1/2, 1/3, 1/5 and 
1/10. Quantitative differences between dynamical protocols 
are seen even at the highest temperature. Data correspond 
to the mean of 10 trajectories; error bars are shown sparsely 
for clarity. 



as they assemble. We use the virtual-move algorithm (en- 
forcing a diffusion constant scaling D(n) oc n^ 1 ) to calcu- 
late the likelihood of the correlated motion of a subset of 
size n of a given cluster, normalized by the likelihood of 
monomer unbinding from that cluster. At high tempera- 
ture single-particle unbinding is rapid on the timescale of 
whole-cluster motion, while at low temperature the con- 
verse is true. However, even at the highest temperature 
we find that motions correlated on timescales less then 
At contribute appreciably to the 'spectrum' of dynamic 
relaxation. Such motion is ignored within single-particle 
algorithms. The non-monotonicity of P move (n) reflects 
the fact that moves of sub-clusters of size n ~ n ma _ x /2 
involve the breaking of many bonds, and so are amongst 
the least favorable processes in an energetic sense; the 
asymmetry of P move (n) results from the fact that clus- 
ter diffusion constants decrease with increasing n. Each 
data set in Figure [5] was obtained by applying 10 7 moves 
to a cluster of size n max (recording any accepted move 
but not making that move) , and averaging over 10 3 such 
clusters. 

In Figure [6] we show lattice gas configurations at fixed 
time following a quench from high temperature to a low 
final temperature, for four values of this final temper- 
ature and for two particle concentrations (here and in 
subsequent sections we used VMD [37] to render simula- 
tion configurations). We employ both collective dynam- 
ics a = 1 and single-particle dynamics a = 00. The equi- 
librium states ('equil.') are in all cases phase-separated 
configurations, which we deduce by using nonlocal Monte 
Carlo moves [35]. At the highest temperatures, coarsen- 
ing dynamics are visually similar. In this regime the sys- 
tem lies 'outside' the clustodal for all but the smallest val- 
ues of n, and single-particle evaporation and diffusion is 
the dominant mode of relaxation. At lower temperatures 
and the lower of the two particle concentrations, strik- 



ing differences emerge: for sufficiently low temperatures 
collective modes of motion (a = 1) induce kinetic frustra- 
tion through the mutual collisions of clusters that bind 
awkwardly and fail to relax before encountering similar 
structures. If collective motion is suppressed (a = 00) 
then instead isolated, compact structures are formed. 
At the higher particle concentration, steric effects par- 
tially frustrate large-scale collective motion, and coars- 
ening patterns are qualitatively more similar (although 
differences can be clearly seen). 

We quantify the influence of collective motion upon as- 
sembly dynamics by measuring the characteristic length- 
scale of domains as a function of time at <fio = 0.1, for 
a = 00 and a = 1. We calculate the domain lengthscale 
I from the first moment of the structure factor [55]; we 
display in Figure [7] results for four temperatures. At the 
highest temperature, T — 1/2, single-particle and cluster 
algorithms show the same qualitative behavior, with both 
I and dljdt increasing monotonically with time. Both 
dynamics produce compact clusters, as may be seen in 
Figure [6] However, even at this high temperature we 
observe quantitative differences between the two algo- 
rithms, with collective motion giving rise to larger do- 
mains at any given time than does single-particle motion. 
The mean cluster size does not become large enough that 
whole-cluster diffusion is negated as a viable means of as- 
sembly. Interestingly, and counter-intuitively, the growth 
in lenghscale under the two algorithms is more similar at 
late times at a slightly lower temperature of T = 1/3. 
Here we observe the same qualitative assembly behavior 
under the two algorithms as at T = 1/2, but now the 
mean cluster size is larger at late times because of the 
stronger thermodynamic driving force. Even when col- 
lective dynamics is permitted, collisions involving these 
larger clusters are sufficiently rare that increases in £ are 
driven chiefly by single-particle unbinding events; the two 
algorithms behave similarly in this regime. 

At the two lowest temperatures, T — 1/5 and 1/10, 
single-particle dynamics displays the same qualitative be- 
havior as it does at T = 1/2 and 1/3, but cluster moves 
show different behavior. While I still increases mono- 
tonically with time, the slope dljdt has at late times a 
constant value at T — 1/5, and decreases at T = 1/10. 
The corresponding pictures for (j>o = 0.1 in Figure [6] sug- 
gest that this decrease in slope results from the steric 
hinderance associated with percolating fractal-like clus- 
ters. Such clusters coarsen chiefly by way of rare single- 
particle unbinding events, a mechanism that leads to 
smaller increases in domain length with time than do 
mutual cluster-cluster collisions. 

Our results show that the 2d attractive lattice gas dis- 
plays a rich variety of coarsening or self-assembly behav- 
iors, and that these behaviors depend sensitively upon 
the dynamical protocol used to evolve the system. Fig- 
ure [6] demonstrates that when one accounts for collec- 
tive motion (a = 1), nucleation and growth mechanisms 
(T = l/2,0o = 0.1) yield to gelation as temperature 
is reduced (T = 1/5 and 1/10, <fio — 0.1). Gelation at 
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low temperature is driven by explicit collective motion, 
or motion correlated on the fundamental timescale At: 
when we forbid such motion (a = oo), gelation is avoided 
in favor of very slow nucleation and growth. It is inter- 
esting that collective motion appears to be crucial to the 
formation of large-lcngthscalc gels within the attractive 
lattice gas at low particle concentration, but that such ex- 
plicit correlated motion is not required to observe glassy 
behavior within a model of (dense) silica [2U]. These 
results imply a difference in the nature of the dynami- 
cal cooperativity associated with glasses and attractive 
gels, respectively, short-time cooperativity versus long- 
time cooperativity. 



we also compare 'freely-draining' with 'Stokesian' clus- 
ter diffusion scalings. Freely-draining motion implies 
(D trans oc n~ ,-Drot. where n is the number 

of monomers comprising the cluster and / is the cluster 
moment of inertia about the rotation axis. For Stokesian 
scalings we take (-Dtrans. oc R~ 1 ,D T at. oc i?~ 3 ), with R a 
measure of the cluster radius of gyration perpendicular 
to the translation vector or axis of rotation. 



- A Brownian 

- A Stokes 
■A SPM 

■ B Brownian 

- B Stokes 
-B SPM 



IV. THE ATTRACTIVE LATTICE GAS 
GENERALIZED TO CONTINUOUS SPACE: 
MODEL ASSOCIATING COLLOIDS 

The effects of collective motion seen in the lattice gas 
are also apparent in a simple off-lattice model of hard 
spheres equipped with isotropic interactions. Such a 
model is a caricature of colloidal particles in solution 
with small polymers; the colloids associate by virtue of 
a polymer- mediated depletion attraction. In Ref. 39J, 
attractive colloidal particles of this nature were observed 
in experiment to assemble into clusters whose geometry 
depended upon the range and strength of the depletion 
attraction. Strikingly, at low colloid concentration these 
clusters appeared to be stable, defying the expectation 
that components with strong pairwise attractions should 
phase-separate or gelate. 

We consider as a simple model of this system a three- 
dimensional collection of N hard spheres of diameter a 
equipped with an attractive square well of range £cr and 
strength U. Spheres are placed randomly within the sim- 
ulation box and occupy 4% of its volume. We impose 
periodic boundary conditions. We consider two of the 
parameter sets described in Ref. [39] : a potential of mod- 
erate strength and range, A (U — 2.6 k^T, £ = 0.11), and 
a potential of considerable strength and short range, B 
(U = 12fc B T,£ = 0.02). We carried out Monte Carlo 
simulations of systems of N = 1085 and 6500 particles, 
drawing particle displacement magnitudes (in units of 
a) uniformly from the interval [0,0.15], and rotation an- 
gles uniformly from a distribution with maximum ~ 14° . 
We used 'virtual- move' translations (employing the al- 
gorithm described in the appendix), and used a static 
linking scheme to effect rotations about the center of 
mass of a chosen pseudocluster. These choices of dis- 
placement magnitude and rotation angle imply a basic 
timescale At sufficiently large that to a good approxima- 
tion large-scale cooperative motion cannot occur from 
uncorrelated moves of single particles; such moves ef- 
fect only local structural relaxation and the binding, un- 
binding and diffusion of monomers. Our aim is to com- 
pare such dynamics with pathways accessible to explic- 
itly correlated motion. When we consider cluster moves 
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FIG. 8: Kinetic data for the self-assembly of model col- 
loids. We show number of clusters n as a function of Monte 
Carlo time t for systems A (U = 2.6 k B T, £ = 0.11) and B 
(U = 12fc B T,C = 0.02), under 'freely-draining' (Brownian) 
and 'Stokesian' cluster diffusion scalings. System A assem- 
bles or coarsens chiefly through binding, unbinding and diffu- 
sion of monomers, and displays little dependence upon cluster 
diffusivity. These data superpose on data generated using a 
single-particle Monte Carlo algorithm (SPM). System B in- 
stead assembles through the concerted motion of aggregates 
(when collective motion is permitted), and displays a depen- 
dence upon cluster diffusivity. When denied collective motion, 
clusters coarsen via very slow monomer unbinding events. 
Data are averaged over 5 stochastic trajectories, each gen- 
erated using 1058 particles. Error bars are displayed sparsely 
for clarity. 



As in the lattice gas, collective motion plays a dom- 
inant role when potential energy gradients encountered 
on a timescale At are large, and is less important than 
single-particle unbinding events for small energy gradi- 
ents. In Figures [8] and [9] we show kinetic measures for 
systems of 1058 particles. System A displays no statis- 
tically significant difference in the evolution of the num- 
ber of clusters with time if collective motion is allowed 
or suppressed. System B, however, experiences dramat- 
ically different fates under collective and individual par- 
ticle motion, with the former inducing gelation and the 
latter resulting in a slowly-coarsening collection of iso- 
lated clusters. 



In Figures 10 and 11 we display snapshots of these 
systems for collections of 6500 particles (in these images 
system A was evolved using a maximum translation of 
0.3<t). System A undergoes phase separation into crys- 
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FIG. 9: Kinetic data for the self-assembly of model colloids. 
We show the number of clusters of given set of sizes Si, n Si , as 
a function of Monte Carlo time t for system B under Brownian 
and Stokesian cluster diffusion scalings, and for single-particle 
moves (SPM). Here S\ — 3, S2 denotes cluster sizes from 7 to 
10, and S3 denotes cluster sizes from 26 to 64. Dynamical 
pathways for the two realizations of collective motion differ 
even for relatively small aggregates. Data are averaged over 
5 stochastic trajectories, each generated using 1058 particles. 
Error bars are displayed sparsely for clarity. 



talline clusters. The effect of collective motion is appar- 
ent only at late times when clusters fuse; such fusing is in 
qualitative agreement with Brownian dynamics simula- 
tions of a similar system [40 . System B forms under col- 
lective motion stringy, frustrated aggregates that merge 
and form a gel, with the timescale for gelation different 
for the two different cluster diffusivities modeled. In ex- 
periment, colloids with a strong depletion attraction (of 
nature similar to systems A and B considered here) form 
instead isolated clusters and do not phase-separate 
or gel. The authors of Ref. [40 discuss possible reasons 
for the disparity between this experimental observation 
and the cluster-cluster aggregation seen in simulations in 
Ref. [40] (similar to those observed here): one such sug- 
gestion is that, in experiment, an accumulation of charge 
on large bodies might induce a repulsion that stabilizes 
a phase of isolated clusters. 



V. A MODEL OF VIRAL CAPSID 
SELF-ASSEMBLY 

The ability of interacting proteins to spontaneously 
form icosahedral capsids in vivo and under some con- 
ditions in vitro is a striking example of biological self- 
assembly. Protein subunits assemble under a variety 
of conditions, avoiding both kinetic and thermodynamic 
traps. Understanding the mechanisms that render as- 
sembly so robust is an essential step towards designing 
synthetic analogs of viral capsids. In addition, a com- 




FIG. 10: Structures obtained from Monte Carlo simulation 
of 6500 particles equipped with attractive interactions of type 
A (top row, excluded-volume view) and B (bottom row, bond 
view), permitting (left column) or forbidding (right column) 
explicit collective (Stokesian) motion. The behavior of system 
A differs under the two dynamical protocols only at late times, 
where collective motion leads to the aggregation of crystalline 
clusters. For system B, collective motion induces gelation 
while single-particle motion results in small, slowly-ripening 
clusters. These behaviors are similar to those of the lattice 
gas (see Figure IHJ. Times of image capture (clockwise from 
top left), in units of 10 6 MC steps, are 1.4, 2.9, 1.7 and 0.6. 




FIG. 11: Assembly pathway for system A under explicit 
collective moves. Assembly is driven principally by single- 
particle binding, unbinding and diffusion, with the effects of 
collective motion apparent only at late times when phase- 
separated clusters fuse. Times of image capture, in units of 
10 6 MC steps, are from left to right 0.2, 0.6 and 1.4. 



prehensive understanding of the viral capsid formation 
process will spur the design of antiviral drugs and new 
drug delivery systems; the latter could possess the ability 
to assemble and disassemble around their cargo without 
requiring explicit external control. 

In Ref. [9] a class of simple models was introduced 
in order to study the mechanism by which interacting 
protein subunits might form 'capsids', 60-member closed 
shells having icosahedral symmetry. We focus here on 
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the 'B3' model of that reference. Units interact via a 
pairwise potential that models an excluded volume and 
a short-ranged, angularly specific attractive interaction 
designed to stabilize capsids. 

Here we examine the role of collective motion within 
the assembly dynamics of this model. We evolved 
via virtual-move Monte Carlo a collection of 1000 sub- 
units. These are initially randomly oriented and dis- 
persed within a three-dimensional simulation box; we 
take periodic boundaries in each dimension. We scaled 
cluster diffusivities in order to approximate Brownian dy- 
namics. As found in Reference [9], we observe a regime of 
inter-unit potential strength and specificity within which 
assembly is robust. We show in Figure [l2]the capsid yield 
(fraction of units residing in complete capsids) at fixed 
observation time as a function of potential strength eb 
or specificity (9 m . Yield data obtained via virtual-move 
Monte Carlo agree with Brownian dynamics results to 
within statistical error (considerable variations in yield 
are observed within each algorithm at a given thermo- 
dynamic state). Particle displacement magnitudes are 
drawn from a uniform distribution with maximum equal 
to a length unit er; particle interaction range is 2. 5c Ro- 
tations are scaled accordingly. Yields are non-monotonic 
functions of these parameters, for the reasons outlined in 
Figure [T] overly strong or insufficiently specific interac- 
tions promote malformed intermediates that fail to as- 
semble into complete structures; overly weak or specific 
interactions result in productive subunit-subunit bind- 
ing events that are too rare to induce assembly on the 
timescales simulated. We show in Figure fl3] example con- 
figurations obtained from a well-assembled and a badly- 
assembled system. 

We can examine the role of collective motion in this ex- 
ample of self-assembly by explicitly restricting or forbid- 
ding collective modes of relaxation within Monte Carlo 
dynamics. We find that in many cases final yields are not 
strongly affected by doing so, but that assembly dynam- 
ics are impaired. In Figure [12] we show yields obtained 
using single-particle Monte Carlo dynamics (triangles) at 
Monte Carlo times equal to (dotted lines) the times at 
which collective-move Monte Carlo data were sampled, 
and for times at which yields appeared to saturate (dot- 
dashed lines). Saturated yields were not obtained (after 
300 hours of simulation) for single-particle moves in the 
right panels of Figure [12] 



In Figure 14 we present measures of assembly kinetics 
obtained at two parameter sets for fully collective motion 
(we move clusters according to the algorithm described 
in Section ITT] with the diffusivity of n-mers chosen to be 
-D(n) oc n in order to approximate Brownian dynam- 
ics), for motion allowing explicit moves of monomers and 
dimers only (as for fully collective motion, but with the 
diffusion constant of trimers and higher-order clusters set 
equal to zero) and for single-particle motion (diffusion 
constant of dimers and higher-order structures set equal 
to zero). To generate these data we drew particle dis- 
placement magnitudes from a uniform distribution with 




FIG. 12: Final yield y of capsid model as a function of 
potential specificity # m (left column) or strength £b (right 
column). Left panels: yield at subunit concentration Co = 
Na 3 /L 3 = 0.11 for attraction strength 16 k B T (top left) or 
14 k B T (bottom left); Right panels: yield at binding speci- 
ficity # m = 0.5 for concentration Co = 0.037 (top right) or 
Co = 0.11 (bottom right). Circles denote Brownian dynam- 
ics results [9]; squares denote results obtained using collective 
Monte Carlo dynamics; triangles denote results obtained us- 
ing single-particle Monte Carlo dynamics at times equal to 
(dashed lines, right column) those at which collective-move 
data were sampled, or at later times (if available) when yield 
had appeared to saturate (dot-dashed lines, left column). 
Times of data capture for Monte Carlo simulations: left pan- 
els, 4.2xl0 6 MC steps; right panels 9xl0 6 MC steps. Excep- 
tions are the late-time single-particle data captured at 9xl0 6 
MC steps (top left panel, triangles) and 7.1 x 10 6 MC steps 
(bottom left panel, triangles). 




FIG. 13: Configurations of the capsid model generated using 
virtual-move Monte Carlo. Left: configuration illustrating 
high yield obtained at parameter set eb = 16A:bT, 6> m = 0.5, 
Co = 0.11. Right: malformed shell of 76 particles obtained at 
parameter set eb = 22 feT, # m = 0.5, Co = 0.037. 



maximum equal to 0.2c. While at these thermodynamic 
states the final yields do not depend strongly upon the 
availability of correlated moves, the kinetics of assembly 
is more rapid when collective motion is allowed. Interest- 
ingly, restoring explicit moves of dimers alone is sufficient 
to recover much of the ease of assembly afforded by a 'full' 
spectum of correlated motions. These results support the 
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FIG. 14: Capsid kinetics at binding energy eb = 14 k^T and 
concentration Co = 0.11 evolved using Monte Carlo dynam- 
ics with different degrees of collective motion permitted. We 
show yield y and energy per particle e as a function of time. 
Assembly under collective moves is more efficient than under 
single-particle moves. Restoring collective motion of dimers 
recovers a substantial fraction of the efficiency of fully collec- 
tive motion. Data in left panel are the mean of 10 trajectories, 
data in right panel are the mean of 4 trajectories. Error bars 
are displayed sparsely for clarity. 

observations made in Ref. [9]. There it was found that 
B3 capsids grow in part through events consisting of col- 
lisions between intermediates larger than monomers, im- 
plying that collective motion plays an important role in 
the model's assembly dynamics (see also Ref. [Hj). Fur- 
ther, the most frequent intermediate binding events for 
the B3 model involved dimers. Indeed, we find here that 
explicit dimer motion renders assembly much more facile 
than if moves are uncorrelated. We note also that the 
tendency of subunits to form closed shells suppresses in 
large part the kinetic trapping seen in the other models 
studied in this paper, where large aggregates bind in an 
awkward fashion and frustrate equilibration. 

These capsid assembly results may be contrasted with 
the assembly properties of a model of interacting pro- 
tein complexes called chaperonins, studied in Ref. |30j . 
There the non-complementarity of model chaperonin- 
chaperonin interactions coupled with collective modes of 
motion induce a high degree of kinetic frustration in some 
regions of the phase diagram: large clusters bind awk- 
wardly, which slows or prevents equilibration. Single- 
particle moves, which suppress such cluster diffusion, give 
rise to small, isolated, well-formed structures. We con- 
clude that particle interaction geometry and the extent 
of organized aggregates play a decisive role in shaping the 
effects of collective motion in self-assembling systems. 

VI. CONCLUSIONS 

We have examined the role of collective motion in 
Monte Carlo simulations of three model systems. We 
find that collective motion is responsible for gelation at 
low particle concentrations within the two dimensional 
lattice gas, at thermodynamic states for which the equi- 



librium configurations are phase-separated. Suppressing 
collective motion at the same thermodynamic states in- 
stead leads to nucleation and growth of clusters. Col- 
lective motion plays a similar equilibration-frustrating 
role within a simple off-lattice model of associating col- 
loids in three dimensions, driving gelation when inter- 
actions are strong. By contrast, correlated motions of 
anisotropically- interacting subunits within a model of vi- 
ral capsid assembly lead to more efficient self-assembly 
under all conditions considered. Our results suggest that 
gelation of homogeneously-interacting particles might be 
regarded as coarsening in the presence of collective mo- 
tion, and that the interplay of phase-separation and gela- 
tion depends both upon the thermodynamic state, and 
upon the rate of diffusion of self-assembled aggregates. 

Control of collective motion in real systems is very dif- 
ficult, but might be possible in special circumstances. In- 
dividual colloids in high-salt polymer solutions interact 
via only short-ranged depletion attractions, but buildup 
of charge on large aggregates can lead to many-body elec- 
trostatic repulsions. It is possible that fine-tuning of solu- 
tion conditions could be used to select the smallest aggre- 
gate lengthscale for which repulsions become significant, 
thereby controlling the extent to which multimers collide 
and bind. Many-body effects of a different nature can 
be induced by long-ranged hydrodynamic interactions in, 
for example, scdimcnting colloidal suspensions [42] . In 
the case of viral capsids, productive multimer-multimcr 
bindings occur when aggregates' exposed contacts meet; 
changing the number of contacts per subunit (by muta- 
tion, for instance), thereby changing the sticky surface 
presented by a capsid of a given size, might be used to 
affect the importance of multimer binding to assembly 
pathways. 
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VIII. APPENDIX A 



choice of pij (/i — > f) , 



The virtual-move scheme discussed in Ref [SU] and the 
main text accounts for the asymmetry of link formation 
at the level of the acceptance rate. Here we present a 
modification of this scheme in which link formation is 
forbidden if the corresponding link would not form un- 
der the reverse move. We find that the acceptance rate 
for collective motion is made simpler. Consider first 
a generic dynamic pseudocluster-formation procedure, 
moving from state /i to state v, in which a link is formed 
between particles i and j with probability Pij(fi — > v). 
This probability is computed by making a virtual move 
of i. We subsequently specialize to the particular choice 
of (/j, — > v) given in Equation ([2| . 

1. Start in state /j,. Choose a seed particle, say i, and 
a move map. Add % to the pseudoclustcr, the list 
of particles to be moved. 

2. Choose a neighbor of i not in the pseudocluster, 
and with which i has not in the current move pro- 
posed a link. Call this neighbor j. With probability 
PijifJ, — > v) (computed by moving i under its vir- 
tual map) form a 'pre-link' between i and j (not an 
automatic link, as described in the main text). 

• If the pre-link does not form, we consider that 
a link has failed to form (outright failure). 
Choose another neighbor of i, say k, and re- 
turn to stage 2, with the replacement j — > k. 



If the pre-link ij forms, compute 
reverse link balance factor /reverse 
mm 1, t - iJA — 



the 



1 Pijin^u) 

With probability /reverse, form a full link 
between i and j. Add j to the pseudo- 
cluster. Go to particle j and proceed from 
step 2, with the replacement j — > j' and 
i -> j. 

If a full link fails to form, record the link 
ij as frustrated. Do not add j to the pseu- 
docluster. Choose another neighbor of i, 
say k, and return to stage 2, with the re- 
placement j — > k. 



3. Proceed until no more links remain to be tested, 
and evaluate the acceptance probability for the 
move. 

The acceptance probability for the revised algorithm fol- 
lows from a modification of Equation (JTJ in the main text. 
We choose to balance the total rates for forward and re- 
verse moves involving a given realization 1Z of internal 
pseudocluster links, and a realization of internal failed 
links 'blind' to the nature of those failed links (whether 
frustrated or outright failed) . The acceptance rate for the 
dynamic linking procedure described here, for a generic 



W^ U) = @{n c -n c )V(C) 
x min|l,e" /3( - E "" B " ) 

YCt^ mi^ -> m) qui" -> m) 

X 

x K1B) 



(ij) 



PijiH^u) min(l, ^{^ 



The first two lines of Equation ( 15 1 are as Equation ([IJ 
of the main text. Variables q denote outright failed 
links between C and its environment. When we take 
PijifJ- v) as in Equation ^ of the main text, such 
variables cancel the Boltzmann bond weights for all but 
a specialized class of moves (see main text). Variables 



Pij 



denote 



qij{n -+v)= pij{fi -> v) (l 

frustrated links between C and its environment. Frus- 
trated links cannot form between the pseudocluster and 
its environment for both forward and reverse moves, and 
should such links form during the forward move then that 
move must be rejected. Variables q denote unformed 
links internal to the pseudocluster, whether frustrated 
or outright rejected. We have that 



qij(fi ->i/) = 



qij(fi — ► v) + qij{(i -> u) 
1-Pij(fj, -> v) 



Pij(fi -> v) I 1 - mi 
1 — Pij([i — > v) min 



1 Pij( v 



M) 



Pij(n -»■ v 
1 Pijj v -> M) 



,(16) 



and so ?ij(/i v)jq\j[y — > u)=l. Lastly, the product in 
the final line of Equation ( |15[ ) runs over all fully-formed 
links. By construction of the linking procedure each quo- 
tient in this product is unity: we ensure that links formed 
during the forward move can also form during the reverse 
move. 

If we take Pij(fA — * v) as in Equation (2| of the main 
text, than the acceptance rate for the collective move is 



W&*"W = Q(n c -n c )V(C) 



x mm 



hSi.c. n 



.(.17) 



The link-formation and link-failure factors internal to the 
pseudocluster have canceled. The factor 6{. B . in the third 
line of Equation ( 17 1 is unity if no frustrated links join 



the pseudocluster to its environment, and zero otherwise. 
This factor is required because a frustrated link indicates 
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that a move of one particle relative to another has effec- 
tively been rejected. If we then do not form the link, in 
order to see if both particles are incorporated into the 
pseudocluster, we must reject any move in which both 
particles do not end up in the pseudocluster. The ad- 



vantage of this scheme relative to that presented in the 
main text is that here links internal to the pseudocluster 
that form under the forward move, but do so with zero 
probability under the reverse move, do not automatically 
result in the rejection of that move. 
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